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ABSTRACT 

The  Sharp-Wheeler  bubble  model  of  Rayleigh-Taylor  instability  is  stu- 
died through  numerical  computations.  By  using  as  input  the  proper  coefficient 
of  the  single  bubble  velocity  and  the  coefficient  of  the  bubble  merger  cri- 
terion, and  by  setting  the  initial  bubble  distribution  based  on  the  dispersion 
relation  of  linear  Rayleigh-Taylor  instability,  we  find  a  numerical  solution  for 
the  bubble  interface  as  a  function  of  time.  The  interface  is  observed  to  be  in 
almost  constant  acceleration  with  an  acceleration  coefficient  a  approximately 
equal  to  0.059.  The  results  are  in  qualitative  agreement  at  least  with  the 
experiment  of  Read. 

1.    Introduction 

We  describe  in  this  paper  a  detailed  analysis  of  Rayleigh-Taylor  instability  in  its  non- 
linear and  chaotic  stage.  In  this  stage,  the  Rayleigh-Taylor  instability  is  characterized  by  the 
motion  of  bubbles  and  spikes.  According  to  a  theory  of  Sharp  and  Wheeler',  rwo  types  of 
motion  will  determine  the  dynamics  of  the  bubble  interface:  the  single  bubble  motion  and  the 
merger  of  two  bubbles. 

There  are  two  dimensionless  constants,  c,  and  c;,  which  characterize  these  motions  and 
hence  the  dynamics  of  the  Sharp-Wheeler  model.  Solving  this  model  numerically,  we  find 


where  h^,  is  the  mean  position  of  the  bubble  interface.  A  scaling  law  shows  that 


a  =  a-,4-.  (2) 


We  find  that  a^-,  is  a  nearly  universal  constant  in  the  sense  that  it  depends  only  on  the  initial 
bubble   distribution   and   is   nearly   a  constant   for   a  significant  range  of   initial   data.    Here 

A  =  — ^ is  the  Atwood  number  of  the  flow. 

P:~Pi 

There  is  a  body  of  knowledge  concerning  a,  c;  and  C;.  In  this  paper,  we  tie  some  but 
not  all  of  these  facts  together,  using  the  Sharp-Vr  heeler  model. 

Constant  acceleration  (1)  of  the  bubble  interface  was  reported  by  Read-.  We  reproduce 
this  behavior  numerically.  Read  reports  that  a  is  a  universal  constant,  independent  of 
materials,  and  we  show  a  range  of  parameters  for  which  a,-  is  a  universal  constant  in  (2). 
This  shows  a  qualitative  agreement  of  the  Sharp-Wheeler  model  with  Read's  experiment. 
However,  there  are  problems  with  the  coefficients  c,  and  ct- 

First,  c-  has  been  studied  extensively  numerically,  analytically  and  experimentally  for 
single  bubbles.  In  this  body  of  work,  c^/^A  has  a  universal  value  approximately  equal  to 
0.32  in  two  dimensions.  We  find  that  Read's  experiments  do  not  agree  with  this  value.  It 
appears  that  c,  in  the  context  of  a  chaotic  flow  has  to  be  renormalized  from  its  value  for  reg- 
ular, single  bubble  flow.  Moreover  for  chaotic  flow,  cy'^ A  does  not  appear  to  be  a  universal 
constant.  <:2  has  not  been  studied  in  detail.  Our  preliminary  numerical  calculations  of  c;  indi- 
cate an  order  of  magnitude  agreement  with  Read's  experiment. 

We  take  experimental  values  for  c,  and  C;  from  Read,  and  insert  them  into  the  Sharp- 
Wheeler  model.  The  result  is  in  good  quantitative  agreement  with  Read's  measured  value 
a  =  0.062.  Because  of  the  lack  of  understanding  of  ci  for  chaotic  flow,  the  incomplete 
understanding  of  c;  for  any  flow,  and  the  fact  that  in  only  a  single  experiment  of  Read  could 
c;  be  measured,  this  quantitative  agreement  is  provisional.    Furthermore,  we  note  that  the 
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maximum  position  /r^,,  diverges  more  rapidly  from  h-  m  the  Sharp-Wheeier  model  that  it 
does  in  the  experiment,  so  that  quantitative  agreement  cannot  apply  to  all  aspects  ci  the  flow. 

On  the  basis  of  dimensional  analysis,  the  single  babble  has  a  terminal  velocity  which  is 
proportional  to  the  square  root  of  gravity  times  the  bubble  radius.  It  can  be  written  as 

i 
h   =   cigr)-.  (3) 

where  h  is  the  bubble  height,  g  is  the  gravitational  constant  and  r  is  the  bubble  radius.  We 
will  summarize  estimates  of  the  coefficient  r,  by  various  investigators  using  experimental, 
analytical  and  numerical  methods.  We  will  also  present  a  determination  of  c  from  Read's 
experiments. 

The  bubble  merger  process  has  been  observed  in  experiment.  It  is  postulated'  that  adja- 
cent bubbles  merge  when  one  is  sufficiently  far  ahead  of  the  other.  After  the  merger,  there 
is  a  single  bubble,  whose  cross  section  equals  the  union  of  the  cross  section  of  its  two  consti- 
tuents and  whose  height  is  determined  by  conservation  of  mass  within  this  cross  section.  The 
criterion  for  merger  to  occur  in  the  Sharp-Wheeler  model  is  the  inequality 

/i,-/!,   2:  m(r),  (4) 

where  h.,  i  =  1,  2  are  the  heights  of  the  two  bubbles  before  merger,  r;  is  the  radius  of  the 
smaller  bubble  and  m(r)  is  an  unknown  function  of  r.  As  a  simple  approximation  we  sup- 
pose that  m(r)  can  be  written  as 

m{r)  =  c.r.  (5) 

In  section  3,  we  will  analyze  C;  from  the  numerical  simulation  using  the  front  tracking  hydro- 
dynamics code  and  the  experimental  results  of  Read-. 

The  two  parameters  c^  and  c^  obtained  from  the  analysis  of  the  single  bubble  motion 
and  the  two  bubble  merger  process  are  then  used  in  a  statistical  computation  in  which  a  large 
number  (about  10")  of  bubbles  evolve  under  the  rwo  given  types  of  bubble  motion.  The  aver- 
age height  of  bubble  interface  as  a  function  of  time  is  found  from  this  computation.    The 


.  a  . 

average  heigh:  h  of  the  computed  bubble  interface  is  fit  by  a  linear  lunction 

h  =  a  (.V  -  Y.-.) ,  (6) 

where 

X  =  A^r-.  P) 

where  t  is  time,  and  X^  is  a  constant.  Both  X  and  X    have  units  of  length.    This  notation  and 
the  definition  of  X  and  a  agree  with  that  of  Read. 

2.    The  determination  of  c^ 

The  single  bubble  motion  for  Atwood  number  A  =  1  has  been  studied  by  several  inves- 
tigators including  Taylor-'.  Birivhoff-,  and  Garabedian'.  For  a  three  dimensional  bubble,  Tay- 
lor gave  the  value  of  c.  =  OAS.  By  using  a  conformal  mapping  method.  Birkhoff  and  Carter 
find  that  the  rwo  dimensional  bubble  has  an  approximate  analytic  solution  which  gives  the 
estimate  c.  =  0.325.  However,  the  solution  was  shown  non-unique  later  by  Garabedian.  He 
showed  that  different  c_  ma\  be  obtained  for  bubbles  with  different  radii.  He  postulated  that 
the  fastest  bubble  (largest  c ■)  would  be  physically  correct.  The  best  estimated  upper  bound 
for  the  values  of  c,  by  Garabedian  is  c,  s  0.34. 

There  have  also  been  numerical  studies  of  the  velocity  coefficient  <"■  of  the  single  bubble 
motion.  The  vortex  simulation  by  Baker  et  al^  for  incompressible  flow  starts  with  a  small 
amplitude  sine  mode  which  evolves  to  bubble  motion  in  the  later  stage  of  a  Rayleigh-Taylor 
instability.  The  bubble  velocity  approaches  a  constant,  and  the  ratio  of  h,'{gr)-  -  is  shown  to 
be  about  0.32.  Menikoff  and  Zemach  did  a  numerical  study  of  the  bubble  motion  using  a 
conformal  mapping  method.    Their  value  for  c^  is  about  the  same  as  that  by  Baker  et  al. 

By  applying  the  front  tracking  hydrodynamics  code,  Glimm  et  al'  have  studied  the  bub- 
ble motion  of  the  Rayleigh-Taylor  instability  in  the  case  of  compressible  fluids.  The  terminal 
velocity  v,^^^^  varies  with  the  Atwood  number  as  well  as  with  the  non-dimensional  compressi- 
bility 


.V/  .  \/  M£.  (8) 
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where  \  ib  the  characteristic  bubble  radius,  v  is  the  gas  constant,  p  and  p  are  the  density  and 
pressure  of  the  lower  gas  respectively.  In  the  incompressibie  case  (.Vf  =  0),  it  is  expected 
that  c  /'^  .4  IS  a  constant,  given  by  rlie  /A  =  1  value  of  c  ^  A  ~  0.32.  It  was  shown  numeri- 
cally t  at  c-.^'A  is  a  monotone  increasing  function  of  M  and  a  possible  A  dependence  for 
\f  s  i  IS  not  resolved  clearly.  The  results  are  shown  m  Fig.l.  According  to  Fig.l,  Ci'^.A 
varies  between  0.32  and  0.63  with  different  Atwood  numbers  and  compressibility.  For 
M  —  0.  the  results  agree  with  the  incompressible  results  cited  above. 

Experimental  work  to  verify  r.  has  also  been  reported  by  several  authors.  Taylor,  in 
supporting  his  theoretical  result,  showed  that  the  3D  bubbie  has  a  constant  ve!ocit\  with 
c,  =  0.49.  The  early  experiment  by  D.  Lewis^  is  somewhat  controversial.  Although  he 
observed  the  constant  velocity  of  the  bubbles,  his  experimental  value  of  c,  is  substantially 
larger  than  the  theoretically  p'edicted  value.  According  to  him,  for  a  bubble  in  water,  c;  is 
about  1.11,  and  an  air-glycerine  surface  has  an  even  larger  c  between  3.5  and  5.5.  In  distinc- 
ton  to  other  authors,  who  use  a  cross  sectional  radius  for  r,  Lewis  used  the  radius  of  curva- 
ture of  the  bubble  in  the  calculation  of  c  .  This  would  account  for  only  a  small  part  of  the 
difference  in  the  reported  values  of  c  . 

The  experiment  done  by  H.  W.  Emmons  et  aP'^  shows  a  much  better  agreement  with 
theory.  In  their  experiment,  an  initially  still  glass  tank  containing  liquid  and  air  is  accelerated 
through  a  vertical  rail  track.  In  their  later  stage,  the  bubbles  are  observed  to  have  a  constant 
velocity.  According  to  Emmons  et  al,  the  measured  c-  in  his  experiment  lies  between  0.28 
and  0.42. 

Ratafia"  used  the  rwo  materials  octyl  alcohol  and  water,  which  give  a  much  smaller 
Atwood  number  of  0.095  (density  ratio  1.22).  The  value  of  Ci/^-A  is  between  0.276  and 
0.457.  This  experiment  stresses  an  important  fact  that  although  c,  varies  among  materials 
with  different  Atwood  numbers,  c.f^A  is  a  universal  constant. 


Al!  of  these  studies  concern  either  a  sineie  bubble.  ;.e.  column  of  light  fluid,  confined 
in  a  tube  or  a  periodic  repetition  of  this  pattern.  In  the  periodic  case,  we  can  insert  iinaginary 
no-flow  boundaries  between  the  individual  bubbles  to  confine  the  flow,  and  thus  regard  the 
flow  as  equivalent  to  the  single  bubble  case. 

In  a  chaotic  flow,  the  speed  of  the  bubble  interface  is  set  by  the  fastest  (largest)  bubble. 
This  bubble  will  move  ahead  of  the  others,  and  for  this  reason  it  moves  with  different  side 
boundary  conditions.  The  leading  bubble  is  not  (or  is  less)  confined  laterally  by  its  neigh- 
bors, and  it  rises  as  a  free  bubble,  unconfined  ;n  a  large  body  of  heavy  fluid.  The  problem  of 
free  rise  of  a  single  bubble  of  limited  extent  has  also  been  studied.  The  corresponding  value 
of  c,  is  2,'3  m  three  dimensions  (Taylor  and  Davies  -,  also  .AUred  and  Blount'M,  which  is 
50"^  greater  than  the  corresponding  value  for  regular  smgie  bubble  flow.  On  this  basis,  we 
see  that  there  is  no  requirement  for  the  idealized  c  of  single  bubble  flow  to  equal  that  of 
chaotic  flow.  We  also  note  that  the  single  free  bubble  is  unsupported  by  the  light  fluid 
below.  Thus  its  c  also  does  not  need  to  coincide  with  the  c,  of  chaotic  flow. 


We  present  our  own  analysis  of  r.  based  on  the  experiment  by  Read.  Five  cases  in 
Read's  experiment  are  selected  to  calculate  cr.  three  of  them  are  in  two  dimensions  and  two 
are  in  three  dimensions  We  find  values  of  c  ./^  A  significantly  larger  than  the  above 
numbers,  for  regular  single  bubble  flow. 

Read  provided  /i  as  a  function  of  X .  which  is  given  by 


X  =  A 


|(.(0)'^ 


dt' 


(9) 


Taking  the  square  root  of  (9),  and  then  it,  derivative  with  respect  to  time,  we  obtain 


X  =  2  [AXg]  ' 


(10) 


The  experimental  data  of  Read  indicates  a  linear  depencjence  of  h  on  X, 


h  =  a(X-.Yo), 


(11) 


and  the  constants  a  and  X^  can  he  found  from  the  experimental  data.  Therefore  .'i  =  aX 
=  2a\AXg\  -  .  Also,  according  to  our  assumption,  h  =  c,g'-r'-.  Comparing  these  two 
formula,  we  have 


AX 


(12) 


The  values  of  r  and  X  are  obtained  from  Read  in  the  following  way.  The  frame  of  the 
pictures  in  Read  are  150  mm  x  150  mm.  We  take  the  radius  of  the  largest  bubble  in  each  pic- 
ture as  the  effective  bubble  radius  for  the  calculation  of  (12;,  excluding  the  bubbles  at  the 
edges  of  the  tank,  which  appear  to  be  exceptional.  This  convention  was  also  followed  by 
Read  in  the  experimental  measurement  of  h.  For  each  picture,  h  is  measured  and  the 
corresponding  X  is  then  found  from  the  same  paper  of  Read.  The  following  tables  show  our 
results  for  c,  as  well  as  c,/^  A  . 


Experiment  34 


a   =   0.062  A   =  0.997 


(2D) 


t(ms) 

r 

h 

.Y 

^1 

r,/^'7  ) 

23.6 

6.35 

15.25 

230.4 

0.74 

1 

0.74     ; 

1 

1 

29.2 

9.79 

21.60 

332.0 

0.72 

0.72 

34.8 

13.23 

30.50 

460.0 

0.73 

0.^3 

-  ^  - 


Experiment  36 


Experiment  29 


1     a   = 

0.064 

.4   = 

0.5 

2D1 

1 

<x   = 

0.0-3 

A    = 

0.997 

(- 

5D)            j 

I 

i 

1  t(ms') 

'.       r 

h 

X 

"' 

|c„\7 

1 

!  t(ms) 

1 

r 

^ 

i       -^ 

'' 

c,/Vl  • 

\r.6 

3.43 

11.44 

232.0 

0.74 

1      1.05 

1 

I    30.8 

i 

3.06 

15.63 

'    244.6 

1.30 

1.30 

i 

1   44.1 

i 

;   5.72 

i 

21.10 

383.0 

0.74 

■      1.05 

1 
1 

t 
1 

;    35.8 

1 

1 

!   4.33 

25.35 

I    314,9 

1.24 

1.24 

1 

1   60.8 

:   6.99 

i 

41.94 

709.0 

0.91 

1 

1 
1 

40.8 

[ 

!   6.33 

t 
! 

33. SO 

i   402.9 

i 

1.16 

1.16 

Experiment  56 


Experiment  35 
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In  summary,  we  conclude  that  aiihouEh  the  velocity  of  the  single  bubble  motion  has 
been  given  by  an  analvtic  theor\',  by  numerical  simulation  ana  experimen..  this  value  of  r. 
does  not  fit  the  data  trom  Read's  experiments.  The  disagreement  can  be  understood  by  the 
fact  that  theory  and  earlier  experiments  involve  oniy  one  or  few  bubbles,  while  the  interface 
in  Read's  experiment  con'ains  many  hubbies.  Therefore  we  speculate  that  the  chaotic  nature 
oi  the  flow  in  Read's  experiment,  whicn  distinguishes  it  from  other  earlier  experiments  and 
theory  could  be  a  cause  of  this  discrepency.  In  fact,  the  calculation  of  c-  from  experiment  37 
for  the  3D  bubble  with  initially  imposed  perturbation  shows  a  better  agreement  with  theory, 
because  in  this  case  the  bubble  is  periodic.  The  c  measured  from  this  experiment  is  about 
0.50  to  0.55  comparing  to  theoretical  value  of  0.-16  b\'  Taylor. 

3.    The  determination  of  r- 

It  has  been  postulated  that  two  adjacent  bubbles  will  merge  and  finally  form  a  single 
bubble  with  a  radius  equal  to  the  sum  of  the  original  two  bubbles  and  the  height  determined 
by  the  conservation  of  mass  within  its  cross  section.  The  merger  is  governed  by  the  criterioa 
(4)  which  means  that  as  soon  as  the  large  bubble  is  a  distance  ahead  of  the  smaller  one  by 
m{r^),  the  two  bubbles  will  merge.  Equation  (5)  can  be  regarded  as  the  first  term  in  the 
Taylor  expansion  for  m(r'). 

No  analytic  theory  has  been  given  to  describe  the  two  bubble  merger  problem,  due  to 
the  complexity  of  the  problem.  Therefore,  C;  has  to  be  determined  by  experiment  and  by 
numerical  study.  The  front  tracking  hydrodynamics  code  is  again  used  to  simulate  the  prob- 
lem. The  initial  interface  configuration  with  two  bubbles  (one  large  and  one  small)  is  created 
by  the  superposition  of  Fourier  modes  in  the  input  data  of  the  computation. 

Fig.  2a  shows  the  successive  motion  of  two  bubble  merging.  The  compressibility  and  the 
Atwood  number  in  this  run  are  M-  =  0.2  and  A  =  0.67  respectively.  .As  we  can  see  from  the 
plots,  the  bubble  with  the  larger  radius  moves  faster  than  the  smaller  one.  .Merger  occurs 
when  the  difference  in  height  of  the  bubbles  is  about  twice  the  radius  of  the  smaller  bubble. 
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The  merger  process  continues  for  a  period  of  time.  According  to  these  pict'^res.  it  is  reason- 
able to  estimate  that  the  coefficient  r-  of  the  merger  criterion  satisfies  c-  ==  2.  Fig.  2b  and 
Fig,  2c  are  the  density  and  pressure  contour  plots  of  the  bubble  merger  respectively.  Fig.  2d 
shows  a  distinct  merger  process  for  the  two  cases  .'►/-  =  0.1,  A  =  0.43  and  M-  =0.1  and 
A  =  0.90  respectively.  Again,  we  have  c-  =  2. 

We  can  also  estimate  c;  from  an  e.xperiment  b\  Read.  Fig.  3  in  Read-  shows  the  inter- 
face of  the  two  fluids  at  two  time  instants.  According  to  these  pictures,  one  bubble  disap- 
peared at  ;  =  34.gmj.  while  the  radius  of  the  other  bubble  increased  by  1/2  of  its  original 
value,  indicating  that  a  merger  occured  between  r  =  29.2m?  and  t  =  34.8m5.  By  measuring 
the  difference  in  heights  of  the  two  bubbles  and  the  radius  of  the  smaller  bubble  at 
f  =  29.2mi,  we  estimate  that  r-  ~  2.  Refined  numerical  studies  of  c-  are  needed  to  complete 
the  discussion  given  here. 

4.    Similarity  and  scaling  law  of  the  statistical  bubble  motion 

Since  at  the  initial  time  of  the  Rayleigh-Taylor  instability,  the  bubble  heights  and  radii 
are  randomly  distributed,  the  overall  motion  of  the  average  bubble  interface  has  to  be 
described  statistically.  The  Sharp-Wheeler  mode!  combines  two  major  features  of  the  bubble 
motion.  I.e.  the  single  bubble  motion  with  constant  velocity  and  the  two  bubble  merger  under 
the  criterion  (4).  We  use  the  code  originally  written  by  McBryan  et  al  -  to  simulate  this  sta- 
tistical model. 

There  is  a  geometric  similarity  in  the  statistical  motion  in  the  Sharp-Wheeler  bubble 
model.  This  can  be  seen  by  reducing  Eq.  (3),  (4)  and  (5)  to  a  set  of  dimensionless  equations. 
We  introduce  rwo  length  scales  in  the  bubble  model.  These  are  the  horizontal  length  scale  R 
which  we  set  equal  to  the  initial  average  bubble  radius  R  =  r,, .  and  the  vertical  length  scale 

H  =  c;/?.    We   nondimensionalize   Eq.    (3),   (4)   and   (5)   by   making  r    =  -^,   /i,  =   -j^,   and 

t  =   -ip,  where  the  characteristic  time  scale  is  7"  =  — -— -.    In  this  wav,  Eq.  (3)  can  be 

T  c^igRy- 


rewritten  as 

4  =  -■   -  (13) 

at 

and  the  merger  criterion  f4)  and  (5)  become 

>i;    -    h-    2    r    .  (14) 

Eq.  (13)  and  fH)  are  dimensionless  equations  with  no  free  parameters  that  govern  the  bub- 
ble evolution.  The  average  radius  and  height  of  the  initial  bubble  distrib'jtion  are 
r_.  =  h..  =  1.  The  only  remaining  variable  aspect  of  the  problem  is  its  initial  distribution 
profile.  In  our  simulation,  the  bubbles  at  the  initial  time  are  uniformly  distributed  between 
r--r  ,  and  r^  -r.,  as  well  as  between  h.  -h  and  n.  -  h  _,  .  We  characterized  the  initial 
bubble  distribution  by  the  two  dimensionless  parameters  f  ,.  =  r^^.'R  and  /;..  =  h...  H.  In 
doing  this,  we  ignore  the  dependence  of  the  distribution  on  its  higher  moments. 

N\'e  can  easily  derive  a  scaling  law  from  the  nondimensionalized  bubble  equations.  The 
equations  tell  us  that  the  parameter  a  introduced  by  Read  is  proportional  to  cf,',4c^.  This  can 
be  proved  as  follows: 

«  =   -^  =   ^^^   =    #-a„  (15) 

For  the  remainder  of  this  section,  we  use  the  reduced  variables  and  omit  the  hat  on  f  etc. 
This  is  equivalent  to  taking  c,  =  C;  =  1. 

An  important  feature  of  the  dimensionless  bubble  model  shown  by  the  numerical  simu- 
lations is  that  the  dimensionless  ratio 

«o  =  -;r  (16) 

is  almost  a  constant  of  0.24  as  long  as  r_.,.  is  not  too  small.  In  Eq.  (16),  h^  is  the  height  of 
the  average  bubble  interface  at  the  time  t.  This  is  true  statistically  in  the  sense  that  the  more 
bubbles  considered  initially,  the  more  nearly  constant  oq  is. 


The  remaining  problem  is  the  dependency  of  a,  on  r ,,  and  h  ..  .  The  numerical  simula- 
tion shows  that  the  change  of  a  '.s  smaii  with  tne  change  of  h_^.  as  we  can  see  from  Fig.  3. 
To  show  that  the  computation  of  a,,  is  statistically  correct,  we  have  also  changed  the  bubble 
ensemble  by  inserting  distinct  seeds  in  the  random  number  generator  and  by  making  further 
random  shu  fling  of  the  random  number  generated  by  the  computer.  The  results  are  all 
presented  in  Fig.  3,  which  show  no  significant  change  in  the  a.,.  We  have  also  shown  that  olq 
varies  little  with  the  change  of  r.,  as  long  as  r..  >  0.25  (Fig.  4).  However,  for  r,,  <  0.25, 
the  linear  dependency  of  h...  on  :-  is  seriously  affected  as  we  can  see  from  Fig.  5.  In  fact,  the 
much  smaller  r...  case  is  similar  to  the  situation  in  which  the  interface  initially  has  a  periodic 
(sine)  wave  perturbation  so  that  all  bubbles  in  the  nonlinear  stage  have  an  almost  equal 
radius.    See  Glimm  et  al''  where  this  phenomenon  was  also  observed. 

We  have  also  tested  several  special  cases  of  the  bubble  model.  The  first  special  case  is 
one  for  which  the  initial  bubble  number  density  has  an  exponential  decay  with  the  bubble 
radius.  In  such  a  case,  we  have  observed  that,  although  the  linear  dependence  of  a^  on  r-  still 
remains  valid,  the  computed  a^  is  substantially  larger  than  in  the  case  of  the  bubble  distribu- 
tion of  finite  radius.  This  result  is  thought  to  be  due  to  the  existence  of  the  few  larger  bub- 
bles in  the  tail  of  the  initial  bubble  distribution  which  appear  to  make  the  merger  process  fas- 
ter. 

The  second  special  case  is  the  two  mode  distribution.  For  this  case,  we  specify  rwo 
separate  groups  of  bubbles  in  the  initial  bubble  distribution.  The  first  group  of  bubbles  has  a 
smaller  r,,  and  relatively  larger  r^^.  and  hj,. .  The  total  number  of  bubbles  in  this  group  is 
much  larger  than  the  second  one  (about  100  :  1).  This  ensemble  represents  the  naturally  gen- 
erated bubbles  in  the  later  stage  of  the  Rayleigh-Taylor  instability.  The  second  group  of  bub- 
bles has  larger  r^,  and  much  smaller  deviation  of  r  and  h.  This  ensemble  represents  an  artifi- 
cially imposed  perturbation  on  the  initial  bubble  interface. 

It  appears  that  for  this  bimodal  distribution,  the  linear  dependence  of  ac  on  t-  no  longer 
exists.    Fig.    6   shows    the    result    of    the    numerical   simulation.    According    to    Fig.    6,    the 


-  13  - 

acceleration  rate  of  the  bubble  interface  at  the  bei»:nning  of  the  evolution  is  much  larger  than 
the  single  mode  result.  This  phenomenon  can  be  characterized  as  a  runaway  situation,  in 
which  the  larger  bubbles  in  the  second  ensemble  quickly  swallow  the  smaller  bubbles  in  t;ieir 
neighborhood.  As  the  small  bubbles  gradually  disappear  and  the  smaller  bubbie  ensemble 
shrinks,  the  acceleration  rate  a  resumes  its  normal  rate.  Fig.  7a  and  Fig.  "b  show  the  com- 
parison of  the  bubble  merger  pattern  between  the  single  mode  distribution  and  the  two  mode 
distribution. 

5.    Comparison  with  experiments  of  Read 

In  order  to  simulate  the  experiment  by  Read,  we  have  first  to  estimate  the  initial  distri- 
bution of  bubble  radii  and  heights.  The  estimation  is  based  on  the  dispersion  relation  of  the 
linear  theory  of  Rayleigh-Taylor  instability.  The  growth  rate  is  thought  of  as  specifying  the 
probability  for  bubble  occurrence.  The  number  density  of  bubbles  is  assumed  to  be  uni- 
formly distributed  about  the  average  bubble  radius  r^.  with  the  deviation  width  r  ,^  ,  where 
r^.,  is  of  the  same  order  of  magnitude  as  the  wavelength  of  the  most  unstable  linear  mode  and 
tj,.  is  estimated  from  the  width  of  the  unstable  modes  of  the  dispersion  relation,  h^^  is 
specified  in  a  similar  manner. 

The  dispersion  relation  is  from  S.  Chandrasekhar'"  and  the  formula  was  first  derived  by 
Harrison  in  1908.  The  dispersion  relation  has  more  recently  been  studied  by  Menikoff  et 
al'-.  In  their  study,  Menikoff  et  al  used  a  nondimensionalized  wave  number,  growth  rate  and 
surface  tension.  For  a  given  surface  tension,  there  are  two  curves  which  describe  the  upper 
and  the  lower  bounds  of  the  solution  of  the  dimensionless  dispersion  relation.  The  variations 
of  the  exact  solution  due  to  the  change  of  other  parameters  such  as  densities  and  viscosities 
are  found  within  11%  by  these  upper  and  lower  bounds.  Due  to  its  complexity,  the  exact 
solution  has  to  be  found  numerically.  For  this  reason,  we  have  written  a  numerical  program 
which  solves  the  dispersion  relation  through  a  bisectional  root  finder.  Fig.  8  shows  the  exact 
solution     with     its     upper     and     lower     bounds.      The    effective     gravity    in    this    case     is 
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g  =  3.98X  l0-(cm/5ec-). 


whicn  IS  equivalent  to  the  effective  gravity  m  experiment  3-1  of  Read. 

The  physical  value  of  the  maximum  growth  rate  at  k  =  20.1(cm~ ■)  indicates  that  the 
most  probable  or  the  average  bubble  radius  initially  is  r  ,  =  0.16cm.  We  can  also  estimate 
the  width  of  the  distribution  of  bubble  radii.  According  to  the  dispersion  relation,  we  have, 
approximately,  r.,.   =  0.01cm. 

With  the  input  parameters  c,  =  .7  and  c^  =  2,  we  have  simulated  the  experiment  34  of 
Read.  Fig.  9  shows  the  bubble  interface  motion  as  a  result  of  the  statistical  bubble  code.  The 
middle  line  plots  the  height  of  the  average  bubble  interface  vs.  gt-  while  the  upper  and  lower 
lines  plot  the  the  heights  of  bubbles  of  maximum  and  minimum  radii  vs.  gt-.  The  horizontal 
axis  unit  of  gt-  is  equivalent  to  X/A  in  Read's  paper. 

The  graph  shows  both  qualitative  and  quantitative  agreement  with  the  experiment  by 
Read.  First,  the  linear  dependence  of  average  height  of  the  bubble  interface  on  X  has  been 
proved.  Physically,  it  means  that  the  average  bubble  interface  is  constantly  accelerating. 
Second,  the  numerically  calculated  acceleration  coefficient  a  is  about  0.059.  which  agrees  with 
the  experimental  value  0.062  of  Read,  well  within  the  uncertainty  in  the  values  for  c,  and  c^. 
The  ratio  /j^.^'/i^    does  not  agree  with  the  experiments. 

As  we  have  shown  in  the  previous  section,  the  constant  acceleration  rate  uq  breaks 
down  when  the  deviation  of  initial  bubble  distribution  T^,.lr^.  <  0.25.  Therefore  there  is  a 
question  as  to  whether  a  so  called  self-resonating  material  interface  exists  in  which  the  bubble 
radii  of  the  Rayleigh-Taylor  instability  initially  have  an  almost  equal  size.  We  can  show  that 
this  is  impossible.  According  to  the  dispersion  relation  of  the  linear  theory,  for  invicid  fluids 
with  finite  surface  tension,  the  cut-off  wavelength  for  the  unstable  modes  is 


\.  = 


^(P:  -  Pi) 


1 : 


(17) 


where  cr   is  the  surface   tension,   while  the  wavelength  of  the  most  unstable  mode   (which 
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represents  the  average  bubble  radius]  is  given  by  ,\.,    =  ^  3\  .    Therefor.,  the  ratio  r^^ir^ 
can  be  estimated  as  toilows 

I-L  ^  A>L  =    ^'-^'    =   1  -   ^  =  0.42.  (IS) 

r.  \...  A..  i 

The  inclusion  of  viscosity  will  not  vary  the  result  much  as  we  can  draw  from  the  upper  and 
lower  bound  theory  by  Menikoff  et  a!'.  This  conclusion  is  a  partial  explanation  as  to  why 
the  acceleration  rate  for  Read's  experiments  is  a  universal  constant. 

6.    A  discussion  of  the  bubble  merger  process 

Eq.(4),  which  describes  merger  as  an  instantaneous  process,  is  an  oversimplified  pic- 
ture. There  has  been  no  systematic  study  of  the  two  bubble  interaction  in  Rayleigh-Taylor 
instability,  neither  theoretically  nor  experimentalh .  The  numerical  simulation  of  the  bubble 
merger  process  by  the  authors  is  still  in  progress.  However,  there  are  several  points  which 
deserve  discussion  in  this  paper. 

First,  the  conservation  of  mass  m  the  Sharp-Wheeler  merger  model  could  be  ques- 
tioned. The  model  states  that  the  radius  of  the  bubble  after  merger  equals  to  the  sum  of  the 
the  radii  of  the  two  bubbles  before  merger  and  the  height  yields  a  conservation  of  mass. 
This  gives  a  counter  intuitive  reverse  flow  of  the  bubble  height.  However,  it  is  observed  from 
the  numerical  simulation  that  the  assumption  of  mass  conservation  is  not  strictly  valid.  Fluid 
above  the  smaller  bubble  is  entrained  into  the  mixed  region  below  the  bubble  interface  dur- 
ing merger  and  this  entramment  contradicts  strict  conservation  of  mass  above  the  bubble 
interface  during  merger. 

Secondly,  the  numerical  simulation  shows  a  rather  complicated  bubble  structure  after 
the  merger.  It  has  been  observed  that  the  merger  is  accompanied  by  a  pinch-off  of  the  larger 
bubble  due  to  the  interaction  of  the  larger  bubble  and  spikes  resulting  in  a  closed  bubble  and 
the  multilayer  structure  of  the  fluids.  To  model  this  multilayer  motion  after  the  bubble 
merger,  which   is   similar   to   the   slug   flow   regime   in   the   multiphase   flow   in   pipes,  could 


require  a  more  complex  model  or  at  least  phenomenological  renormalizar;on  of  c,.  We  also 
note  that  this  pinch-olf  could  be  more  pronounced  in  a  strictly  two  dimensional  computation 
than  It  might  be  in  an  approximately  two  dimensional  experiment. 

Thirdly,  we  note  the  possible  dependence  of  the  bubble  merger  on  such  parameters  as 
the  compressibility  and  the  Atwood  number.  See  also  Fig.  2d. 

Fc  :rthly  there  is  aiso  a  time  dependence  of  the  merger  process.  The  Sharp-Wheeler 
bubble  modei  assumes  that  the  bubble  merger  :s  an  instantaneous  process.  However,  our 
numerical  simulation  shows  that  the  process  takes  a  substantia!  time  period  so  that  the  single 
bubble  motion  and  the  bubble  merger  really  should  be  considered  simultaneously. 

7.    Summary  and  conclusions 

The  bubble  motion  of  Rayleigh-Taylor  instability  is  characterized'  by  rwo  types  of  bub- 
ble movements:  the  single  bubble  motion  with  constant  velocity  and  merger  between  adjacent 
bubbles.  The  value  of  c-  given  by  various  authors  cluster  about  .32  for  two  dimensional  bub- 
bles. The  value  for  C;  appears  to  be  about  2  from  both  numerical  simulation  and  experi- 
ment. 

The  statistical  bubble  code  used  to  simulate  the  evolution  of  a  large  number  of  bubbles 
is  based  on  the  Sharp-Wheeler  model  of  bubble  motions.  It  is  found  that  under  the  merger 
process,  the  average  bubble  interface  is  accelerating.  The  acceleration  rate  a  is  proportional 
to  cr/Ac;.  The  dependence  of  a  on  such  dimensionless  parameters  as  r.^.lr--,  and  h^^Jc-j^.^ 
is  very  weak. 

By  setting  the  initial  bubble  distribution  based  on  linear  theory  of  Rayleigh-Taylor  ins- 
tability, the  solution  of  the  statistical  bubble  code  shows  a  linear  dependence  between  the 
average  height  of  bubble  interface  and  X,  indicating  the  interface  is  accelerated  because  of 
increasing  bubble  radii  through  mergers.  For  parameters  c,  =  0.7  and  c;  =  2,  the  accelera- 
tion coefficient  a  is  equal  to  0.059,  which  agrees  with  the  experimental  value  of  Read  well 
within    the    uncertainty    in    the    value    for    c,    and    c^-    Whether    this    agreement    would    be 


-  17  - 

maintained  ever  a  larger  number  oi  experiments  is  not  known. 

Several  problems  remain.  One  is  to  find  a  :heorc:icai  understanding  of  the  renormaiized 
value  of  C;  appropriate  to  the  chaotic  environment  of  Read  s  experiment.  The  second  is  the 
theoretical  validity  of  the  bubble  merger  model,  which  requires  further  study. 
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Fig.  1  A  plot  of  the  constant  C;  versus  the  square  root  of  the  At^ood  ratio  for 
D  =  1.5,  2  and  3,  \f-  =  .2  (V),  .5  (-),  2  (x  ).  and  y  =  1.4.  The  values  of  C;  calculated 
from  the  t^o  (•  and  A)  dimensional  incompressible  {M-  =  0,  A  =  1)  theories  are  also 
shown.  .At  least  for  small  values  of  At-,  it  appears  that  C;/V^  is  a  constant  depending 
on  .V/-.  The  values  for  M-  =  2  ( x  )  do  not  have  sufficient  precision  to  assert  whether  or 
not  such  a  relation  holds  in  this  case.  This  graph  is  shown  by  the  permission  of  Glimm 
et  al*. 
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Fig.  2a  The  successive  interface  position  of  the  two  bubble  merger  process  from  the  numerical 
results  of  the  hydrodynamic  iront  tracking  code  with  A  =  0.67  and  M-  =  0.2.  The  plots  show 
that  the  large  bubble  on  the  left  merges  with  the  small  bubble  on  the  right.  The  merger  involves 
the  pinch  off  of  the  larger  bubble  and  a  multilayer  fluid  configuration  after  the  merger,  a 
phenomenon  which  has  not  been  considered  in  the  Sheirp-Wheeler  bubble  model. 
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Fig.  2b.  Density  contour  plots  of  the  bubbles  before  and  after  the  merger.  Note  that 
after  the  merger,  there  is  a  multilayer  fluid  confiuration  which  may  affect  the  velocity  of 
the  lager  bubble  after  merger. 
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Fig.  2c.  The  pressure  contour  plots  of  the  bubbles  before  and  after  the  merger.  It  is 
observed  that  during  the  merger,  the  pressure  inside  the  smaller  bubble  is  higher  than 
that  inside  the  larger  one,  a  fact  which  causes  the  splices  to  fall  toward  the  larger  bub- 
ble. 
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Fig.  2d.  Simulations  of  bubble  merger  with  different  Arwood  numbers.  The  case  (a)  is 
the  simulation  with  .4  =  0.43,  and  case  (b)  with  A  =  0.90.  The  compressibihty  Af-  in 
both  cases  is  0.1. 
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Fig.  3.  The  dimensionless  acceleration  rate  a^  vs.  the  normalized  height  deviation  of 
bubble  distribution  at  the  initial  time.  The  numerical  results  show  that  a-,  varies  very  lit- 
tle with  the  change  of  Aj,,  =  >ije,/c;rj, .  The  f,^.  in  this  plot  is  0.5.  Here  the  symbols  O. 
n  and  A  represent  distinct  seeds  for  the  random  number  generator  which  generates  the 
initial  bubble  ensemble. 
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Fig.  4.  The  dimensionless  acceleration  rate  a-  vs.  the  normalized  initial  bubble  radius 
distribution.  The  V  and  A  are  the  approximate  upper  and  lower  bounds  tor  the  varia- 
tion of  Up  in  each  run.  The  results  show  that  a,-)  varies  little  with  the  variation  of  r.,  as 
long  as  r_,,.  >  0.25  is  satisfied.  When  r\,  <  0.25,  the  linear  dependence  of  a,]  on  r 
fails.  The  /i_,.,  in  this  plot  is  0.5. 
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Fig.  5.  The  simulation  result  of  average  height  of  the  bubble  interface  vs.  /-  with  a 
very  small  normalized  deviation  of  bubble  radius  distribution  at  the  initial  time.  In  this 
case,  Tj^,  =  0.1,  and  the  linear  dependence  of  h  on  f-  fails.  The  situation  is  similar  to 
that  of  a  superimposed  sine  wave  perturbation  on  the  bubble  interface  at  the  initial 
time,  so  that  the  bubble  radii  in  their  nonlinear  stage  have  almost  equal  sizes. 
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Fig.  6.  The  height  of  the  bubble  interface  as  a  function  of  r  in  the  bimodal  initial  bub- 
ble distribution.  The  plot  shows  that  the  acceleration  rate  a^  is  much  larger  than  the 
normal  case.  This  is  understood  as  the  runaway  situation.  In  the  early  time  of  of  bub- 
ble evolution,  the  a^  is  about  10  times  of  the  normal  rate.  It  then  gradually  resumes  its 
normal  value  as  the  small  bubble  ensemble  shrinks  and  disappears. 
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Fig.  7.  The  comparison  of  the  merger  pattern  between  single  mode  (Fig.  7a)  and  t\vo 
mode  (Fig.  7b)  bubble  distribution.  The  horizontal  axis  represents  the  x  coordinate  and 
the  vertical  axis  is  the  height.  Each  horizontal  bar  gives  the  space  location  of  the  merger 
event. 
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Fig.  8.  The  normalized  dispersion  relation  with  its  upper  and  lower  bounds.  See 
Menikoff--.  The  physical  parameters  are  selected  to  simulate  the  experiment  34  of 
Read-.  These         are         ^  =  3.95x  iO-t-m/.ifc-,         p;  =  0.0013f/cm,         p;  =  0.78^  cm, 

V-  =  i.S3xlQ'-poises  and  v;  =  \lO~-poises  for  the  middle  curve.  The  lower  bound  is  real- 
ized for  pi  =  p-_  and  v;  =  v;,  and  the  upper  bound  is  realized  for  p;  =  0.    and  v;  =  0. 
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Fig.  9.  A  simulation  of  the  experiment  34  of  Read  by  applying  the  statistical  bubble 
model.  With  c,  =  0.7,  and  c.  =  2,  the  result  gives  a  =  h^.Jgt-  =  0.059,  which  quantita- 
tively agrees  with  Read's  experimental  result  of  o  =  0.062. 
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